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Abstract. To meet the challenge of future deep space programs, an accurate and efficient engineering code for 
analyzing the shielding requirements against high-energy galactic heavy radiation is needed. To address this need, a 
new Green’s function code capable of simulating high charge and energy ions with either laboratory or space boundary 
conditions is currently under development. The computational model consists of combinations of physical perturbation 
expansions based on the scales of atomic interaction, multiple scattering, and nuclear reactive processes with use of the 
Neumann-asymptotic expansions with non-perturbative corrections. The code contains energy loss due to straggling, 
nuclear attenuation, nuclear fragmentation with energy dispersion and downshifts. Previous reports show that the new 
code accurately models the transport of ion beams through a single slab of material. Current research efforts are 
focused on enabling the code to handle multiple layers of material and the present paper reports on progress made 
towards that end. 
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INTRODUCTION 

With its renewed interest in manned space exploration, NASA envisions a number of long duration space missions 
on which astronauts risk exposure to heavy-ion cosmic radiation. Since it has long been recognized that this type of 
radiation can have a detrimental effect on the health of humans and on the performance of mission critical 
equipment, shielding from galactic heavy ions becomes a problem of ever-increasing importance as the space 
program proceeds into an era of extended manned space operations (Space Studies Board, 1996). To meet the 
challenge of future space programs, an accurate and efficient engineering code for analyzing the shielding 
requirements against high-charge and energy (HZE) galactic heavy ions is needed. In the 1980’s, Green’s functions 
techniques were identified as the likely means of generating efficient HZE shielding codes that are capable of being 
validated in laboratory experiments (Wilson et al., 1989). In consequence, a laboratory code, designed to simulate 
the transport of heavy ions through a single layer of material was developed (Wilson et al., 1990; Wilson et al., 
1993). It was based on a Green’s function model as a perturbation series with non-perturbative corrections. The code 
was validated for single layer targets and then extended to handle multi-layer targets (Wilson et al., 1991; Shinn et 
al., 1994; Shinn et al., 1995). This early code used a scale factor to equate range-energy relations of one material 
thickness into an equivalent amount of another material and proceeded to perform calculations of nuclear factors in 
the specific material (Shinn et al., 1997). While this method has proven to be acceptable using low-resolution 
detectors (Shinn et al., 1994; Walker et al., 2005), it is not the most accurate reflection of different material 
properties and is unsuited for high-resolution measurements. Lacking from the prior solutions were range and 
energy straggling and energy downshift and dispersion associated with nuclear events. In a recent publication 
(Tweed et al., 2004), it was shown how these effects could be incorporated into the multiple fragmentation 
perturbation series leading to the development of a new Green’s function code GRNTRN (a GReeN’s function code 
for ion beam TRaNsport). The current version of GRNTRN has proven to be accurate in modeling ion beams for a 



single layer of material (Walker et al., 2005; Tweed et al., 2007), and is now being extended to handle multiple 
layers. Unlike the earlier Green’s function code, the new multi-layer GRNTRN code will not make use of range 
scaling but will instead transport ions through the target layer by layer. 


GREEN’S FUNCTION FOR A SINGLE LAYER TARGET 

The transport of high charge and energy ions through matter is governed by the Linear Boltzmann Equation (Wilson 
1977; Wilson et al. 1991) 

8 x <!> i { X ,E) = Y j \ : a ]k (E,E')^x,E')dE'-a i (E)^ X ,E), X > X \ y =1,2,3,. JV (1) 

k>j 

where N is the number of ions in the model, <j> J ( X ,E) is the fluence of ions of type j moving along the x-axis at 
energy E in units of MeV/amu and cr^E) and <j jk (E,E') are the media macroscopic cross sections. These cross 
sections may be expanded as 

*] (E) = < (E) + a/‘ (E) + a/ (E), a Jk (E, E ') = a/ (E, E ') + a/ (E, E ') + a/ (E, E ') (2) 

where, in each case, the first term refers to collision with atomic electrons, the second term is for elastic nuclear 
scattering, and the third term describes nuclear reactions. Many atomic collisions (~ 1 0 6 ) occur in a centimeter of 

ordinary matter, whereas ~ 10 3 nuclear coulomb elastic collisions occur per centimeter, while nuclear reactions are 
separated by a fraction to many centimeters depending on energy and particle type. This ordering allows flexibility 
in expanding solutions to the Boltzmann equation as a sequence of physical perturbative approximations. 

We are required to solve Equation (1) subject to a boundary condition of the type tfj{ X ',E ) = F (E) . In the case of 
a unit source of mono-energetic k-particles of energy E ' at the boundary, F j (E) takes the special form 

F j (E) = 8 jk 8{E - E '), y = 1,2,3. JV (3) 


and the corresponding solution, which is called the Green’s function, is denoted by the symbol G jk ( X , X ',E,E j . It 

should be noted that the delta function in equation (3) is normalized such that it has units of fluence. Once the 
Green’s function is known, the solution for an arbitrary boundary condition F- -(F) is then given by 

<!>, (x, E) = X \ B G jk (v, X \E,E ") F t (F ") dE " (4) 

k 

In order to construct the solution we rewrite Equation (1) in operator notation by defining the vector array field 
function <t> = \j>. (x, F)] , the drift operator I) = |o ( | and the interaction operator 

i-«— £J. a Jk (E,E')dE'-<jj(E)] with the understanding that I has three parts associated with atomic, 

k 

elastic, and reactive processes as given in Equation (2). Equation (1) is then written as 

[D-r'-r'+CT>o = s r .o (5) 


or equivalently as 


0 = [D-I a, -I e '+CT ] '•Og + J fD-I 
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where is the appropriate boundary condition, Q is the quadrature operator and G° is the solution of the 



homogeneous form of Equation (5) with a unit source at the boundary. G° is called the zero order Green’s function. 
It is associated with the primary fluence and is diagonal matrix operator whose elements take the form (Tweed et al. 
2004) 


GUx,x',E,E') = 


niE'i 


y[27Ts k {x-x\E') Pj[E ] 


exp 


[E-E k (x-x\E')f~ 
2 s k (x-x\Ef 
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where P k (E) = exp 


0 S k (H) 


dH 


is the nuclear attenuation function, S k (E) is the change inf? per unit path 


length per nucleon, R k ( E) is the usual range-energy relation, E k (x -x\E ') = R k | R, (/?')- (x - x ') | is the mean 
energy for incident k-type particles of energy E' after a distance of penetration x — x' (Wilson et al. 1991; pp. 
315), and s k (x - x ', E ') is the corresponding energy width (Wilson et al. 2002). 


Equation (6) is a Volterra integral equation and therefore admits the Neumann series solution 


<D = [G° +Q*G°*E .G° +Q«G°«E 


.Q.G°*E .G° + ...]»® g =[G° +G 1 +G 2 +...]-<t» s 


( 8 ) 


The above formalism lends the following interpretation of the solution. The operator G° propagates the particles 
with attenuation processes. The first term, G° <t> B , propagates the ions at the boundary to the interior. E' -G 0 
is the production density of first generation secondaries at depth x x . These are propagated to the interior by 
G° •E r -G° . Lastly, G 1 4> g = Q G° E' G° ® g represents the sum of all the first generation secondaries 

being propagated from the interval [x',.r] and so on. Since G° has now been identified, the remaining terms in the 
Neumann series (8) may be found via the recurrence formula 

G"=[Q.G°-E r >G" _1 , n> 1 (9) 


Construction of the Neumann series has been discussed in detail in Tweed et al. (2005) where accurate analytical 
approximations are obtained for the first three terms and the remainder of the series is estimated by a non- 
perturbative technique. 


GREEN’S FUNCTION FOR A MULTI-LAYER TARGET 

It will now be shown that the single layer Green’s function described above can be used to generate the Green’s 
function for a multi-layer target. The target is assumed to consist of M layers 

L m={ x:x m- m = 1, 2,3,...,M (10) 


some or all of which may be comprised of distinct, homogeneous materials. The single layer Green’s function for 
the layer L m will be denoted by the symbol G Jk (x,x m _ l ,E,E m ~ l -,rn) and the multi-layer Green’s function by the 

symbol G jt (x,x 0 ,E,E°) . Conceptually, the procedure is quite simple. A mono-energetic beam of k-particles 
fzf (x 0 ,E) = 5 jk 5{E - E° ) enters the first layer, L x , at x 0 , is transported through the layer where fragmentation occurs 
and exits as the multi-particle fluence </> j (x l ,E) = G jk (x l ,x 0 ,E,E°) = G jk (x l ,x 0 ,E,E°;Y) . The fluence exiting Z, 
serves as the boundary condition for L 1 which it enters as if> j (x l ,E) and exits as 


= XI G,(.v...v,/;.//:2)^(.v. //),/// = Xf ^{^^EM^G^^H ,E°)dH (11) 



and so on. Continuing in this way, we see that the fluence at x e L m is given by 


</j(x,E) = Y i l g G JI (x,x m _ i ,E,H-,rn)t i (x m - 1 ,H)dH ^^^x^EM^G^^H ,E°)dH (12) 

i i 

or in operator notation 


<J>{x,E) = G(x,x,„_ l ;»;)*G(x m _ 1 ,x 0 ), x n _, <x<x„, 


(13) 


where, on the right-hand side, the energy arguments have been suppressed in the interest of brevity. It is now clear 
that the multi-layer Green’s function is given by the recurrence relations 


G(x 0 ,* 0 ) = U, G(x,x 0 ) = G(x,x m _ x ;rn)'G(x m _ x ,x 0 ), x m _, <x<x m 


(14) 


where m = 1,2,3 is the layer number and U is the unit or identity operator. Additionally, by expanding each 
term in Equation (14) in a Neumann series, it is readily found that 


't ( 
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J V 9=0 


(15) 


and hence that the n' h order multi-layer Green’s function is given by 

n 

G"(x,x 0 ) = X G (x, x m _, ; m)'G”~ s (x m _, , x 0 ), x m _ t <x< x m 

.s=0 

with n = 0,l,2,3...;w = 1,2,3 


(16) 


The Zero Order Multi-Layer Green’s Function 

Equation (16) shows that the zero order multi-layer Green’s function is given by 

G °(x,x 0 ) = G 0 (x,x m _ 1 ;w)*G 0 (x m _ I ,x 0 ), x m _ x <x<x m 


(17) 


or equivalently by 


G 


: (X, X 0 ,E,E°) = 8 jk G» (X, x,„_, , E, //; m)G° kk (x m _ t ,x 0 ,H,E°)dH, 




(18) 


with m = 1,2,3,...,M. 

The integral in Equation (18) may be evaluated by numerical quadrature or by analytical approximation, the latter 
being preferred in the interest of computational efficiency. In this paper, both methods are used, the numerical 
integration providing a standard of comparison for the analytical techniques developed. 

It is now well known (Tweed et al., 2004) that if a primary beam with a delta function or Gaussian energy profile 
enters a homogeneous layer of material then its energy profile on exit is approximately Gaussian. It may therefore be 
assumed that the primary beam entering the layer L m takes the form 


G^x m _ x , Xo ,E,E°) = 


A m ~\E) 

/--) ~ m - 1 

V2 ns 
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( E-E m ~' V 


(19) 



where the amplitude A'" 1 (E) is a slowly varying function of E , and the mean E m 1 and width s m 1 are independent 
of E . It follows that, for x e L m , 


G° kk (x,x 0 ,E,E°) = j 
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(20) 


The slowly varying terms in this integral will be evaluated at the peak of the second Gaussian and the argument of 
the first Gaussian will be expanded about this peak as follows: 


where 


E k ( x ~ x m-i E k (x - x m _ x , E m ~ l ;m) + r k (x - x m _ x , E"" 1 ; m)(H - E m ~' ) 


r 4 x ~ x m -\ ’ E "' l '- m )=[d H E k (x-x m _ l ,H-,m)] 


H=E m ~ 1 


S k \- E k (x ~ ^,„-i , E '"- 1 ; m ); m] 
S k [E m -'-m\ 


(21) 


(22) 


It can then be seen that 


Gl k {x,x 0 ,E,E°)~ k 


P k [E \m] A m -\E m -') 


-exp^ 


[E-E k (x-x m _ l ,E m -'-,rn) 
s k (x~x 0 ,E ° ) 


P k [E',m] 42ns k (x-x 0 ,E°) 

K. 

where 

s k (x-x 0 ,E 0 ) 2 =s k (x-x m _ i , E"'- 1 ;m) 2 + r k {x - x m _ x , E"'-' ;m) 2 (s"‘-') 2 

It should be noted that the fluence entering L m+l is now given by 


G kk (x m ,x 0 ,E,E °) : 


A"‘(E) 

42ns" 


-exp^ 


f E-E" x2 


(23) 


(24) 


(25) 


where A"'(E) = 


E k \E m ' ; m ] 
P k [E',m] 


A m -'(E m -') , E m = E k (x m -x m _ l ,E m -'-,rn) and s m = s k (x m -x 0 ,E°). 


Since the fluence entering /, (l has amplitude A" = 1 , prescribed mean energy E" and prescribed width .v° . which is 
zero in the case of a delta function input, the above expressions can be used to approximate the primary beam 
fluence in each layer in turn. 

Some results for the zero order Green’s function are presented in Figures 1 and 2. In Figure 1, a comparison is made 
between the single layer theory and the multi-layer theory for a mono-energetic iron beam at 1000 MeV/amu on a 
target consisting of 24 g/cm 2 of Aluminum and a target consisting of three 8 g/cm 2 layers of the same material. The 
figure exhibits graphs of the primary ion fluence at target depths of 0, 8, 16 and 24 g/cm 2 and shows that the multi- 
layer theory is consistent with the single layer theory. 

The target in Figure 2 consists of 8 g/cm 2 of Aluminum followed by 8 g/cm 2 of Polyethylene ( ( CH-, )" ) and 8 g/cm 2 
of Water, respectively, and the beam has a Gaussian profile with mean energy 1000 MeV/amu and initial spread 10 
MeV/amu. The figure presents a comparison between the primary fluence as computed by numerical quadrature and 



as obtained by the analytical approximation described above. The excellent agreement between the results shows 
that the approximate zero order Green’s function accurately depicts the primaiy beam fluence in multi-layer targets 



FIGURE 1. The primary fluence at various depths in 
a 24 g/cm 2 A1 and in three 8 g/cm 2 layers of Al. 



FIGURE 2. A comparison between the numerical 
and approximate primary fluence at various depths in 
8 g/cm 2 Al, 8 g/cm 2 Polyethelyne,8 g/cm 2 H 2 0. 


The First Order Multi-Layer Green’s Function 

By virtue of Equation (16), the first order multi-layer Green’s function is given by the expression 

G) k (x, x 0 ,E,E°) = ^ Gjj ( x, x m _, , E, //; m)G' jk (x m _, ,x 0 ,H,E°)dH + J ? G) k (x, x m _ t , E, H ; m)G^ ,x 0 ,H,E ° ) dH (26) 

Again, the integrals can be evaluated by numerical quadrature or by an approximation procedure similar to that used 
for the zero order case. In this case, the approximation takes the form 

P [E ; m]S ;[E ; m] _ 

G , (v, * 0 , E, E ° ) * 11 ' ' G , (x m _, ,x 0 ,E J ,E°) + A - 1 (E"-' )G , (x, x m _, , E, E”~' ; m) (27) 

Pj[E;m]Sj[E;m] 

where 


E = E (x - x m _, ;m) = R: 1 [i?,. ( E ; m) +x- ; m], x m _ x <x< x m 


(28) 


Figures 3 and 4, respectively, show the first generation i6 0 and 40 Ca fragments from a Gaussian 56 Fe beam of mean 
energy 1000 MeV/amu and width 10 MeV/amu on a target consisting of 8 g/cm 2 of Aluminum followed by 8 g/cm 2 
of Polyethylene and 8 g/cm 2 of Water. Both numerical and analytical techniques are used to compute the fragment 
fluences, which are then exhibited in graphical form for comparison. Again, there is remarkably good agreement, 
which suggests that the analytical approximations made are reasonable. 



The Second Order Multi-Layer Green’s Function 


In this case, Equation (16) yields the expression 


, (x, x 0 ,E,E°) = G°jj (x, x m _ , , E, //; m)G 2 jk (x m _ t ,x 0 ,H,E°)dH + j G] k (x, x m _, , E, //; m)G° kk (x m _, ,x 0 ,H,E°) dH 


+ Z L G J‘ {x ’ ’ E - //; W)G * dH 

i=j + 1 


(29) 


Very good analytic approximations can be obtained for the integrals in the first two terms of Equation (29) but the 
integrals in the last term have not yielded to the techniques used so far. Therefore, for the time being, the second 
order Green’s function G 2 k (x,x 0 ,E,E°) will be approximated by the energy scaling method used by Tweed et al. 
(2004) to construct the non-perturbative Neumann series remainder. To do this, we make use of the integral fluence 

g jm (x-x 0 ) = \ s G Jm (x,x 0 ,E,E°)dE (30) 

which, in a uniform material, is given by the series 

co 

gjmix-x o) = Yjg"m(x-x o) = s jm g(m) + <7 ]m g(j, m) + £ cr jk cr km g(j,k,m) + ■■■ (31) 

n= 0 k 


where the g{j i ,j 2 -,'"ij n ij n+ i-,) are Wilson’s g-functions (Wilson et al., 1989; Wilson et al., 1994). By an argument 
similar to that used in the derivation of Equation (16), it is not difficult to show that the corresponding multi-layer 
quantities g" k (x — x 0 ) are given by the recurrence relations 


gj-k(°)= s jA o’ ^ - *o ) = z (* - *„-i ; «)gr s - x 0 x x m _ t <x< Xm 

i- j 5-0 


(32) 


Assuming that for n > 1 there is little variation in the spectral shapes, we then make the approximation 


G) k {x,x 0 ,E,E 0 )* 8 _* { * G) k (x,x 0 ,E,E°) 

gA x - x o) 


(33) 


The second order Green’s function is illustrated in Figures 5 and 6. 

Figures 5 and 6, respectively, show the second generation 16 0 and 40 Ca fragments from a Gaussian 56 Fe beam of 
mean energy 1000 MeV/amu and width 10 MeV/amu on a target consisting of 8 g/cm 2 of Aluminum followed by 8 
g/cm 2 of Polyethylene and 8 g/cm 2 of Water. Again, the fluence is computed by numerical quadrature and by the 
scaling approximation and then exhibited in graphical form for comparison. 

The scaling method used to approximate the second order Green’s function is clearly not as good as the saddle point 
method used to approximate the zero and first order Green’s functions, but it may nevertheless prove to be sufficient 
for engineering applications. 


CONCLUSION 


A recently developed Green’s function code (GRNTRN) has proved to be very effective in predicting the transport 
of HZE ions through a homogeneous target material. In this paper, steps have been taken to generalize the model 
used by GRNTRN with a view to producing a multi-layer version of the code. Excellent approximations have been 
developed for the first two terms of the Neumann series solution and a less accurate but nevertheless useful 



approximation has been obtained for the third term. Completion of the solution will require the estimation of the 
Neumann series remainder by a non-perturbative technique. Construction of the non-perturbative remainder can be 
based upon the accurate second term approximation, a numerical approximation to the third term or an analytic 
approximation to the third term. All three methods will be developed and the best one used to complete the solution. 
The code will then be validated against ion beam experiments and space based measurements. 



FIGURE 3. A comparison between the numerical 
and approximate first generation of 16 0 fragments in 
8 g/cnr Al, 8 g/cm 2 Polyethylene, 8 g/cm 2 H 2 0. 



FIGURE 5. The second generation 16 0 fragments in 
8 g/cm 2 Al, 8 g/cm 2 Polyethylene, 8 g/cm 2 H 2 0. 



FIGURE 4. A comparison between the numerical 
and approximate first generation of 40 Ca fragments 
in 8 g/cm 2 Al, 8 g/cm 2 Polyethylene, 8 g/cm 2 H 2 0. 



FIGURE 6. The second generation 40 Ca fragments 
in 8 g/cm 2 Al, 8 g/cm 2 Polyethylene, 8 g/cm 2 H 2 0. 
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